#Nature's Kidneys: the Role of Wetland Reserve Easements in Restoring Water Quality
#Nicole Karwowski and Marin Skidmore
#Replication File
#2025

#Link Water Quality Monitor Data to Subwatersheds

###############################################################################
#Set up space

#open libraries
library(fst)
library(ggplot2)
library(raster)
library(sf)
library(dplyr)
library(readr)
library(mapview)


#working directory
setwd("C:/Users/16308/OneDrive - Montana State University/UW_Madison/Water quality/Replication Package")

#open huc 12 data
huc12<-st_read("Data/WatershedBoundaries/Raw/USGS_WBD/USGS_WBD.shp")

#wq sites
sites <- read_csv("Data/SNAPD/Raw/SNAPD/_C_workflow/SNAPD_final_wqd_sites.csv")

sites_sf<-st_as_sf(x=sites, coords=c("x", "y"),
                   crs=st_crs(huc12))

rm(sites)

#make sure they are the same crs
crs(huc12)
crs(sites_sf)

#explore wbd data
summary(huc12)

#cut out unneeded variables
length(is.na(huc12$noncontr00))
length(is.na(huc12$sourcedata))
length(is.na(huc12$sourcefeat))
length(is.na(huc12$sourceorig))
length(is.na(huc12$metasource))
length(is.na(huc12$gnis_id))
length(is.na(huc12$noncontrib))

colnames(huc12)

huc12<-huc12[, c(3,4, 6:9, 11:15, 17, 20)]

table(huc12$hutype)

huc12$areaacres<-as.numeric(huc12$areaacres)
huc12$shape_leng<-as.numeric(huc12$shape_leng)
huc12$shape_area<-as.numeric(huc12$shape_area)
huc12$areasqkm<-as.numeric(huc12$areasqkm)

table(huc12$states)

#do the merge
#use st_within

a<-st_within(sites_sf, huc12, sparse=TRUE)
b<-as.data.frame(a)

sites_sf$row.id<-1:nrow(sites_sf)
huc12$col.id<-1:nrow(huc12)

c<-merge(sites_sf, b, by.x="row.id", by.y="row.id", all.x=TRUE, all.y=FALSE)

d<-subset(as.data.frame(huc12), select=c(col.id, huc12, tohuc, states, areaacres, name))

e<-merge(c, d, by.x="col.id", by.y="col.id", all.x=TRUE, all.y=FALSE)

length(unique(e$huc12))
#there are 100K huc12s in our daat, 21K have site monitors

table(e$states)
colnames(e)

sites_huc12<-e[, c(3:9)]
colnames(sites_huc12)  

#save final data set
save(sites_huc12, file="Data/SNAPD/Processed/sites_huc12.RData")

################################################################################

